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Abstract 


Automating  the  control  of  an  aircraft  flying  in  formation  necessitates  the 
extension  of  the  theory  of  formation  flight  control  to  allow  for  three  dimensional 
maneuvers.  The  formation  was  modeled  as  a  two-aircraft,  leader  and  wingman, 
formation.  Both  aircraft  has  its  own  three  dimensional,  rotating  and  translating,  Cartesian 
axes  system,  with  special  attention  being  given  to  the  motion  of  the  leader  in  relation  to 
the  wingman.  The  controller  operated  using  the  equations  of  motion  expressed  in  the 
rotating  reference  frame  of  the  wing  aircraft.  The  control  system  has  seven  states,  three 
inputs  and  three  disturbance  signals  to  model  the  dynamics  of  the  formation  in  three 
dimensional  space.  The  control  law  employed  was  the  feedback  of  the  difference 
between  in  actual  separation  distance  and  the  commanded  separation  distance  to  affect 
changes  in  thrust,  lift,  and  roll  rate.  The  control  system  incorporated  proportional, 
integral,  and  derivative  control  elements,  each  with  separate  gains,  to  achieve  and 
maintain  the  specified  formation  geometry  despite  various  maneuvers  flown  by  the 
leader.  Simulated  maneuvers  included:  an  initial  displacement  of  the  wingman  away 
from  the  formation  geometry,  and  changes  in  the  leader's  velocity,  altitude,  and  heading. 
For  each  maneuver,  the  controller  performance  was  sufficient  to  maintain  the 
commanded  formation  geometry. 
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THREE  DIMENSIONAL  FORMATION  FLIGHT  CONTROL 


1  Introduction 

Migrating  geese  have  for  centuries  flown  in  formation  to  take  advantage  of  the  air 
flowing  off  of  the  wings  of  preceding  birds  to  reduce  induced  drag.  During  World  War  I, 
when  military  aviation  was  first  taking  to  the  air,  the  United  States  Navy  employed  the 
idea  of  multiple  aircraft  flying  together  for  mutual  benefit  when  it  hung  pursuit  planes 
from  the  dirigibles  AKRON  and  MACON.  The  lessons  learned  during  the  air  campaign 
over  Germany  during  World  War  II,  concerning  the  importance  of  fighter  escort  for 
strategic  bombers,  brought  about  research  programs  such  as  FICON  and  wingtip  coupling 
that  attempted  to  physically  attach  fighters  to  bombers.  Today,  as  the  USAF  turns  to  an 
Air  Expeditionary  Force  mode  of  operations,  the  opportunity  for  mutually  beneficial 
flight  is  ripe  by  applying  the  concept  of  flying  aircraft  in  close  formation.  Pilots  have  the 
ability  to  fly  in  very  precise  formations,  however,  asking  a  pilot  to  accomplish  this  task 
continuously  for  several  hours  is  unrealistic.  Hence,  this  research  seeks  to  further  the 
development  of  an  automatic  close  formation  flight  control  system  for  the  wing  aircraft 
by  extending  previous  work  to  include  full  three  dimensional  aircraft  formation 
dynamics.  An  automatic  formation  controller  which  entails  proportional,  integral,  and 
derivative  control  action  is  designed,  and  its  performance  is  examined  by  way  of  multiple 
simulated  maneuvers. 
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2  Literature  Review 


The  background  references  reviewed  for  this  research  fit  into  three  categories. 
The  first  category  was  an  exploration  of  the  history  of  flying  for  mutual  benefit.  The 
second  pertained  to  aerodynamic  interaction  considerations  for  formation  flight.  The 
final  category  related  to  previous  formation  control  system  investigations. 

2.1  History  of  Flying  for  Mutual  Benefit 

In  reference  [1],  several  concepts  for  flying  together  were  expounded.  The  first 
such  concept  appears  to  have  been  evaluated  by  the  Navy,  who  created  special  attach 
points  on  the  dirigibles  AKRON  and  MACON,  from  which  pursuit  aircraft  were  hung. 
This  combination  of  aircraft  could  loiter  over  hostile  territory  for  reconnaissance 
purposes  and  also  protect  itself  from  enemy  aircraft.  Providing  fighter  escort  for  a 
bomber  flying  long  duration  missions  was  the  impetus  for  the  FICON  (fighter  conveyor) 
project  described  in  [1],  The  conveyor  aircraft  was  the  massive  B-36  outfitted  with  an 
elaborate  trapeze  mechanism  to  which  the  fighter  was  attached  and  then  hoisted  into  the 
bomb  bay.  The  fighter  was  the  XF-85,  an  egg-shaped,  jet-propelled,  folding-wing, 
landing-gearless  experimental  aircraft  that  proved  too  unstable  to  consistently  reattach 
itself  to  the  conveyor  after  flying  its  mission.  The  TOM-TOM  project  followed  up  on  the 
FICON  idea,  but  this  time  the  fighter,  an  F-84,  had  its  wings  remain  outside  of  the  B-36 
fuselage. 
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The  wingtip  coupling  experiments  described  in  [2]  are  interesting  not  only  from 
the  perspective  of  mutually  beneficial  flight,  but  also  from  a  control  system  point  of  view. 
The  principle  underlying  the  wingtip  coupling  concept  was  to  extend  the  range  of  several 
aircraft  by  connecting  them  together  at  the  wingtips,  thereby  effectively  creating  a  single 
aircraft  with  a  very  large  aspect  ratio.  Early  work  with  a  C-47A  and  a  Q-14B  was  very 
successful,  and  the  program  was  extended  to  a  B-29  and  two  F-84  aircraft.  Numerous 
flights  were  accomplished,  but  always  with  the  fighter  pilots  constantly  controlling  the 
flight  of  their  aircraft.  The  next  step  in  the  program  was  an  attempt  to  allow  the  F-84 
aircraft's  automatic  flight  control  system  to  maintain  proper  flight  attitude.  Tragically, 
the  control  system  was  unable  to  properly  adjust  for  the  aerodynamic  effects  of  flying 
coupled  to  the  B-29,  and  both  aircraft  were  destroyed  and  the  aircrews  killed. 

2.2  Close  Formation  Flight 

The  mutually  beneficial  flight  concept  explored  in  this  research  is  called  close 
formation  flight.  The  benefit  comes  in  the  form  of  reducing  the  induced  drag,  achieved 
by  flying  aircraft  in  a  specific  formation  so  as  to  capitalize  on  the  aerodynamic 
interaction  with  the  vortices  emanating  from  the  leading  aircraft's  wings. 

2.2.1  Vortices  Explained 

Vortices  for  aircraft  in  steady,  level  flight  can  be  visualized  as  a  pair  of  horizontal 
rotational  air  flows  that  start  near  the  wingtips.  Vortices  are  defined  by  two  properties, 
vorticity  and  circulation.  Kuethe  and  Chow  [3]  developed  mathematical  expressions  for 
these  properties  in  an  ideal,  inviscid  flow.  Vorticity  is  defined  as  the  rotation  of  a  fluid 
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flow  about  a  point  and  is  measured  in  units  of  length  per  time,  whereas  circulation  is 
defined  as  the  strength  of  the  rotation  and  has  units  of  length  squared  divided  by  time. 
Crow  [4]  proposed  a  stability  theory  for  trailing  vortices,  however,  for  close  formation 
flight,  vortex  instability  is  not  an  issue  because  of  the  relatively  tight  aircraft  spacing.  In 
Kopp  [5],  the  formula  for  velocity  distribution  in  a  viscous  fluid  is  presented.  Rather 
than  having  infinite  vorticity  at  the  axis  of  rotation  as  seen  in  the  inviscid  model,  the 
expression  shows  zero  velocity  at  the  vortex  center.  Figure  1  shows  sample  velocity 
distributions  for  the  two  different  models.  Note  that  as  the  radius  increases,  the  two 
models  yield  similar  answers.  The  velocity  distribution  in  Figure  1  is  found  by  treating 
the  vortex  as  a  filament  at  the  center  of  rotation. 


Figure  1.  Velocity  Distribution  in  Vortex  of  Viscous  and  Inviscid  Flow 
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2.2.2  Vortex  Models 


It  is  known  that  an  aircraft  creates  a  wake  of  vortex  flow  from  each  wing, 
however  some  difficulty  arises  in  modeling  the  flow.  One  such  model  is  called  a  horse 
shoe  vortex,  where  the  wing  is  replaced  with  a  bound  vortex  and  vortex  filaments  extend 
reward  and  are  connected  at  an  infinite  distance  by  another  filament  to  complete  a  circuit 
in  keeping  with  the  Helmholtz  laws.  Other  options  are  to  model  the  aircraft  wake  as  flat 
sheets  of  individual  vortices  or  as  a  rolled  up  sheet  of  vortices  as  demonstrated  by 
Beukenberg  and  Hummel  [6],  Rather  than  an  analytical  model,  Maskew  [7]  developed  a 
numerical  method  called  vortex  lattice  calculations  to  compute  the  flow  properties  of  an 
airplane  wake.  In  this  research,  the  wake  was  modeled  as  a  simple  horse  shoe  vortex  as 
in  Proud  [8]  and  Blake  [10],  allowing  for  the  derivation  of  a  relatively  simple  analytical 
expression  for  the  vortex  flow. 

2.2.3  Formation  Geometry 

The  key  to  reaping  the  benefits  of  close  formation  flight  is  to  fly  in  the  correct 
position  with  respect  to  the  vortex  flow.  There  have  been  several  studies  conducted  to 
determine  the  optimum  aircraft  location  relative  to  one  another  in  the  formation.  The 
optimum  position  for  the  purposes  of  this  report  is  defined  as  producing  maximum 
induced  drag  reduction  while  at  the  same  time  keeping  a  safety  margin  of  longitudinal 
and  lateral  spacing  to  prevent  collisions.  Figure  2  provides  an  explanation  of  spacing 
terminology  essential  to  understanding  the  discussion  of  formation  geometry.  For 
longitudinal  spacing  in  an  echelon  formation,  Maskew  [7]  proposed  three  span  lengths 
was  best  because  that  is  where  the  upwash  reached  a  maximum  value.  Beukenberg  [6] 
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and  Myatt  [11]  both  place  the  lateral  spacing  in  such  a  way  that  wing  tips  actually  overlap 
slightly  to  produce  maximum  upwash  on  the  trailing  aircraft.  The  work  presented  in  [10] 
shows  that  maximum  drag  reduction  is  achieved  when  the  formation  is  all  in  the  same 
vertical  plane  during  steady  level  flight.  Therefore,  in  this  research  the  formation 
geometry  has  a  vertical  spacing  (zc)  of  zero  (referring  to  Figure  2,  zL=0),  a  longitudinal 
separation  distance  (xc)  of  three  times  the  wing  span  of  the  leader  and  a  lateral  separation 
distance  (yc)  of  7t/4  times  the  wing  span  of  the  leader.  It  can  be  seen  that  the  separation  in 
the  x  direction  provides  a  safety  margin  to  counteract  the  wing  tip  overlap  in  the  y 
direction. 


Wing  and  Vertical  Tail  of  W 


Figure  2.  Two  Aircraft  Formation  Schematic. 
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2.3  Automatic  Formation  Flight  Control 

The  work  presented  in  [8],  [9],  and  [12]  modeled  two  aircraft  that  were  fixed  into 
a  single  orientation,  thereby  making  the  controller  design  problem  essentially  two 
dimensional.  The  linearization  of  the  equations  of  motion  allowed  a  proportional-integral 
(PI)  control  system  to  be  designed  using  state  space  techniques.  The  control  law 
developed  in  [8]  was  used  as  the  starting  point  for  the  design  of  the  three  dimensional 
formation  flight  control  system. 
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3  Three  Dimensional  Formation  Dynamics 


3.1  Rotating  Frame  of  Reference 

Formation  flight  control  necessitates  the  use  of  rotating  frames  of  reference. 
Subscript  W  refers  to  the  number  two  or  wingman  aircraft  in  the  formation  and  subscript 
L  refers  to  the  aircraft  leading  the  formation. 

Each  aircraft  is  considered  as  a  point  mass  for  developing  the  dynamic  model. 
Each  aircraft  is  treated  as  a  single  horseshoe  vortex,  with  the  trailing  vortices  separated 
by  7t/4  times  the  span  of  the  wing  and  emanating  from  the  tips  of  the  bound  vortex  that 
simulates  the  wing.  Myatt,  reference  [11],  shows  that  the  lift  distribution  is 
approximately  elliptical,  which  validates  the  use  of  7t/4  times  the  span  of  the  wing  as  the 
distance  between  the  vortex  filaments. 

The  equations  of  motion  are  derived  using  a  translating  and  rotating  frame  of 
reference  at  the  instantaneous  position  of  W.  The  rotating  frame  of  reference  is  a  triad  of 
wind  axes  defined  as  follows:  the  x-axis  aligned  with  the  aircraft  velocity  vector,  the  z- 
axis  aligned  with  the  lift  vector,  and  the  y-axis  (approximately  out  the  left  wing) 
completes  the  right-handed  coordinate  system.  The  derivation  assumes  no  aerodynamic 
side  force  in  the  absence  of  the  L  vortex  filaments,  which  means  no  sideslip  is  allowed,  or 
in  other  words,  coordinated  flight. 
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The  kinematics  are  based  on  the  rotation  from  an  inertial  reference  frame  to  a 
wind  axes  system  attached  to  the  W  or  L  aircraft,  subscripted  according  to  which  aircraft 
is  being  considered. 


Figure  3.  Inertial  and  Wing  Axes. 

The  wind  axes  frame  is  specified  by  the  \j/,  y,  and  <|>  Euler  angles.  For  the  sake  of 
clarity,  the  final  rotation,  about  the  xw  axis  and  through  the  angle  (J),  is  not  explicitly 
shown  in  Figure  3. 

The  orientation  of  the  velocity  vector,  Vw,  of  the  wing  aircraft  is  specified  by 
heading  angle,  \|/w,  and  flight  path  angle,  yw-  The  lift  vector,  Lw,  is  rotated  by  flight  path 
angle,  yw,  and  roll  angle,  <J)W.  An  aerodynamic  side  force  resulting  from  the  vortex 
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filaments  of  L,  Yw,  is  in  the  direction  of  yw,  and  is  included  in  the  derivation  to  allow  for 
the  aerodynamic  interaction  of  close  formation  flight.  Since  a  point  mass  model  is  used, 
moments  are  not  included  in  the  analysis. 

The  three  control  variables  are  lift,  thrust,  and  roll  rate,  and  are  denoted  Lw,  Tw, 
and  pw,  respectively.  The  state  variables  are  Vw  and  the  Euler  angles  \[tw  and  yw- 


3.2  Angular  Rates 

Expressions  for  the  rate  of  change  of  the  three  state  variables  can  be  derived  from 
Newton's  Laws.  Note  that  the  positive  sign  in  front  of  the  weight  term  in  equation  (3.2.1) 
and  the  negative  sign  in  front  of  the  lift  term  of  equation  (3.2.2)  result  from  our  definition 
of  positive  flight  path  angle,  as  shown  in  Figure  3.  Also,  note  that  lift  and  side  force  are 
adjusted  by  bank  angle,  (|>w- 

Vw  =  Tw  ~  +  g  sin  (3.2.1) 

mw 


L\v  cos(j)w  ^gcosy^  Yw  sin  (j)w 
mw^w  tnwVw 


(3.2.2) 


tfw  = - +  Yw  cos<l>w  (3.2.3) 

cosyw  cosy^ 

Resolved  into  the  frame  of  reference  of  the  wing  aircraft,  the  instantaneous 


angular  rate,  a>w ,  of  the  W-wind  axes  triad  is, 


(3.2.4) 


(Ow 


Pw 

Qw 

P\V 


From  [13],  with  flight  path  angle  substituted  for  pitch  angle,  the  above  angular 
rates  are  related  to  the  Euler  angle  time  derivatives  according  to  the  following 
transformation, 


"p" 

'1 

0 

-siny 

< P 

Q 

= 

0 

COS  0 

sin  0  cosy 

Y 

R 

0 

-sin0 

cos 0 cosy 

¥ 

(3.2.5) 


Substituting  the  relationships  for  the  change  in  flight  path 
equation  (3.2.5)  and  completing  the  matrix  multiplication  yields 
for  the  angular  rates  of  the  W-axes  system: 

•  |  tan  y w  (l^  sin  ^  cos  ft, ) 

W  >W  r  r 

mwVw 

Ly,  |  gcosywcos<pw 

y w  ~  T/  +  T/ 
mw  n, 

_  gcosywsm(t>w  Yw 
w  ~  V  mV 

v  w  rnw  v  w 


and  heading  angles  into 
the  following  equations 


(3.2.6) 

(3.2.7) 

(3.2.8) 


3.3  Kinematics 

Fundamental  to  the  problem  of  formation  flight  is  the  relative  position  of  L  with 
respect  to  W.  The  position  and  attitude  of  both  aircraft's  wind  axes  need  to  be  related  to 
an  inertial  reference  frame,  subscripted  I.  Figure  4  shows  the  vector  relationship,  shown 
in  two  dimensions  for  clarity,  between  the  inertial  and  the  rotating  axes  systems. 
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Figure  4.  Rotating  Reference  Frames  (viewed  from  above). 

From  Figure  4,  the  following  vector  expressions  are  developed,  where  the  vectors 
Rl  and  Rw  are  resolved  into  the  inertial  reference  frame  and  the  vector  R  is  resolved  in 
the  rotating  frame  attached  to  W: 

=  ^ 

Therefore  the  rate  of  change  of  the  vectors  can  also  be  written; 

DRl  _  DRW  DR 
Dt  Dt  Dt 

where  the  change  in  relative  position  (R)  in  W's  rotating  reference  frame  can  be 
transferred  to  the  inertial  reference  frame  using  the  following  relationship, 
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(3.3.1) 


DR 

Dt 


dR  _  B 
=  —  +  Q)W  xR 
dt 


The  frame's  angular  velocity,  (dw ,  is  defined  in  equation  (3.2.4).  Therefore  equation 


(3.3.1)  can  be  rewritten  as  follows: 


DRl  _  DR, 


Dt 


Dt 


dR  _  5 

•  +  —  +  (Jt)w  x  R 
dt 


(3.3.2) 


Also,  the  rate  of  change  in  separation  R ,  can  be  written  as  the  time  derivative  of  the 
Cartesian  components  of  the  separation  distance, 


dR 

dt 


xw 

zw 


(3.3.3) 


The  orientation  of  the  rotating  reference  frame  with  respect  to  the  inertial  frame  is 
specified  by  the  \|/,  y,  and  (j)  Euler  angles.  The  three  individual  single-axis-of-rotation 
transformation  matrices  are  in  the  3-2-1  rotation  order  or,  in  other  words,  first  yaw,  then 
pitch,  followed  by  roll. 


cosy/ 

-sinyr 

O' 

cosy 

0 

siny 

1 

0 

0 

*r  = 

siny' 

cosyr 

0 

0 

1 

0 

0 

COS0 

-  sin  (j) 

0 

0 

1 

-siny 

0 

cosy 

0 

sin</> 

COS  0 

The  generic  transformation  matrix,  C,  is  formed  by  multiplying  the  three  single-axis-of- 


rotation  matrices,  viz.. 


C  = 


cos  cos  y 
sin  \j/  cosy 
-siny 


cosy/  sin  y  sin  0  —  sin  yf  cos  (j) 
cos  yr  cos  <j)  +  sin  y/  sin  y  sin  (j) 
cosy  sin  0 


sin  y/  sin  (j)  +  cos  y/  sin  y  cos  0 
-  cos  y/  sin  (j)  +  sin  y/  sin  y  cos  i t > 
cosy  cos 0 


The  preceding  transformation  matrix  is  applicable  to  both  L  and  W,  with  the  only 
difference  being  the  Euler  angles  used.  The  transformation  from  the  L  rotating  reference 


3-6 


frame  to  inertial  is  denoted  as  C[ ,  and  from  the  W  rotating  reference  frame  to  inertial  as 
C'w.  Using  the  transformation  matrix,  C'w ,  the  velocity  relationship  from  equation 
(3.3.2)  can  be  written  in  the  inertial  reference  frame  as. 


cos yL  cosi/AL 

cosy„,  cos  i/a  w 

VL 

cosyLsini/AL 

=  VW 

cos  y  ^  sin  y/  ^ 

-sinyL 

-siny^ 

+  C1 


dR 

dt 


■  +  (0 


\ 

xR 

) 


Solving  the  preceding  equation  for  the  rate  of  change  in  separation  distance  gives, 


dR_ 

dt 


r 

cosyL  cosy/L 

cosyw  cosi/aV(, 

\ 

VL 

cos  yL  sin  i/a l 

-Vw 

cosy„,  sint/A^ 

V 

_  sin  y  L 

-sinyw 

/ 

(3.3.4) 


Using  equations  (3.2.4)  and  (3.3.2),  the  matrix  operations  can  be  performed  to  yield, 


X 

1 

1 

to 

_ i 

cos  yl  cosy/L 

T 

y 

= 

Pyy  Z  P\V  ^ 

+vLcmT 

cos  yl  sin^L 

~VW 

0 

z 

1 

1 

K 

Ol 
_ 1 

-sinyL 

0 

where  the  vector  that  transforms  the  velocity  of  L  into  the  W  reference  frame  is, 


cosyw  cosyL  cost//^  +  sinyH,  sinyL 

cosyL  siny^  sin<pw  cos y/e  +cosyL  cos(()w  smy/e  -cosy^  sinyL  sm(f>w 
cosy^siny^  sin^  cosi/Ae  -  cosy L  cos (j)w  sin  i/a t  -cosyw  sinyL  sm(/)w 


In  the  transformation  vector  above,  the  angle  \[/e  is  the  heading  error,  defined  as 
the  difference  between  the  L  and  W  heading  angles  as  shown  below. 


We  =Wl~Ww 
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3.4  Formation  Model  Equations  of  Motion 

The  pertinent  equations  of  motion  for  a  two-aircraft  formation  model  relate  the 
velocities  and  attitudes  of  the  wingman  and  the  leader  and  the  distance  separating  them. 
Thus,  the  state  variables  in  the  formation  flight  control  system  are:  Vw,  yw,  t|%,  (]>w,  x,  y, 
z,  Vl,  Yl,  and  \jrL  (e  9110),  and  the  number  of  differential  equations  is  ten.  The  three 
control  variables  for  the  controller  are:  Lw,  Tw,  and  pw  (£  913),  where  pw  is  defined  as 
the  time  rate  of  change  of  W's  bank  angle,  <j)w  .  However,  the  creation  of  a  heading  error 

state  as  shown  above  reduces  the  dimension  of  the  system  by  one.  An  additional 
simplification  can  be  made  to  the  control  system  by  declaring  L's  states  exogeneous 
disturbances.  Hence,  the  three  model  disturbances  would  be  L's  velocity,  heading  angle, 
and  the  y  component  of  lift:  VL,  Yl,  and  LLsin<|)L.  In  this  way,  the  original  system  model 
of  ten  states  is  reduced  to  the  following  seven  states:  Vw,  Yw,  <(>w,  t|/e,  x,  y,  and  z  (e  9!7). 
The  seven  differential  equations  of  motion  that  model  the  two-aircraft  system  are  as 
follows: 


T  -D 

_  1W 

V  w 


m 


+  gsinyw 


(3.4.1) 


w 


...  _  Lwcos(j)w-Ywsm(l)w  t  gcosyw 
Yw  ~~ 


mw^w 


K 


(3.4.2) 


w 


K  =Pw 


(3.4.3) 


¥e  = 


sin  (j) ^  ^ 


cosyL  cosy^, 


(3.4.4) 
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(3.4.5) 


XW 


yg  cos  Yw  sin  (j)w  -  zg  cosy^cosf^,  zl^  +  yTv 


+  ■ 


+  VL(cosyw  cosyLcost/^  +  sinyiv  sinyL) 


mw^w 


-v, 


z  tan  y  (l^  sin  0VV  -  Yw  cos  <j)w)-xYw  xg cos yw  cos  0, 
)V  =  3V - + 


mwVw 


V, 


(3.4.6) 


w 


■  VL  (sin  yw  cos  y  L  sin  <j)w  cos  y/e  +  cos  yL  sin  0IV  sin  yre  -  cos  y  sin  yL  cos  (j)w ) 


y  tan  y  (L„,  sin  (pw  -  Yw  cos  (j)w)-xLw  t  xg  cos  y  w  cos  0 

=  -y/v - — j y  +  ^ 


(3.4.7) 


W 


+  VL  (sin  yw  cos  yL  cos  (j)w  cos  y/e  -  cos  yL  sin  (j)w  sin  y/e  -  cos  yw  sin  yL  cos  </>w ) 

The  drag,  Dw,  can  be  eliminated  from  the  velocity  equation  (3.4.1)  according  to 

the  following  derivation.  Using  the  definition  of  the  coefficients  of  lift  and  drag,  an 

expression  for  the  drag  force,  as  a  function  of  lift,  can  be  found. 

D  L 

—  qS 


c 


CL 


where, 


As  in  the  example  on  page  217  of  reference  [14],  the  total  drag  coefficient  can  be  found 
using  the  following  equation: 

CD  =  CDo  +  KC\ 

Substituting  into  the  expression  for  drag  gives. 


D  - 


(  i  \ 


2  A 


+  K\ 


qS 


L_ 

5 s 


3-9 


which  when  applied  to  W  yields, 


Dw  — 


f  r  ^2^ 
\  Qw j 


Qw  S 


w 


Substituting  this  expression  into  equation  (3.4.1)  yields  a  new  velocity  relationship. 


Vw=- 


m 


w 


Tw  qwSwCD 


2  A 


+  gsinyw 


(3.4.8) 


Equations  (3.4.2)  -  (3.4.4)  and  equation  (3.4.8)  entail  the  point  mass  modeling  of 
W  using  wind  axes  coordinate  system  definitions,  and  equations  (3.4.5)  -  (3.4.7)  are  the 
change  in  separation  distances  measured  in  the  W  reference  frame  obtained  from 
equation  (3.3.5). 


3.5  Aerodynamic  Interaction  Effects 

3.5.1  Vortex  Effects  on  the  Wing 

The  effect  of  vortices  on  the  wing  of  W  is  a  function  of  the  separation  distance 
between  the  two  aircraft  and  their  respective  attitudes.  Figure  2  shows  the  two-ship 
aircraft  formation  with  the  vortices  trailing  from  L,  and  the  aerodynamic  surfaces  of  W's 
wing  and  vertical  tail,  all  in  L's  rotating  reference  frame.  The  separation  distance  vector, 

R ,  between  aircraft  can  be  transformed  from  the  W  rotating  reference  frame  into  the  L 
frame  using  the  following  matrix  multiplication: 

Rwl  =  ~C !  Cw  R 


where. 
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c\  =  c{ 


Again  treating  the  wing  of  W  as  a  lifting  line,  as  in  [9],  and  assuming  no  yawing 
of  W,  the  attitude  of  W  is  specified  by  the  y-axis  of  the  W  rotating  reference  frame  and 
can  be  resolved  into  the  L  frame  of  reference.  Points  along  the  wing  line  of  W  can  be 
written  as, 

"O' 

=  -cX  i  i 

o 

Because  W's  wing  is  being  treated  as  a  lifting  line,  the  limits  of  are  the  span 
between  the  locations  of  the  trailing  legs  of  the  horseshoe  vortex,  bvW-  In  the  L  frame  of 
reference,  the  location  of  the  L  generated  vortices  written  in  vector  notation  are: 

V 1  V 

bvL  and  -bvL  Q<rj  < 

°J  [  0 

The  motivation  for  transformation  into  the  L  axes  system  is  this,  the  key  to 
determining  vortex  interaction  is  the  distance  from  the  vortex  to  a  point  on  the  trailing 
wing,  and  that  distance  is  a  function  of  separation  and  attitude. 
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Figure  5.  Separation  and  Attitude  in  the  L  Reference  Frame. 


As  in  [1]  and  [9],  each  wing  can  be  represented  as  a  single  horseshoe  vortex 
system.  Although  the  assumption  of  using  an  elliptically  loaded  wing  to  define  the 
vortex  separation  is  rigorously  true  for  pure  wings  (no  fuselage  effects)  with  high  aspect 
ratio  and  no  sweep  back,  studies  conducted  by  the  Air  Force's  Flight  Dynamics 
Laboratory,  [11],  indicate  that  the  assumption  is  approximately  true  even  for  modem 
fighter  aircraft.  With  this  assumption,  the  effective  separation  of  the  vortex  legs,  called 
bv,  is  7t/4  times  the  wing  span,  b,  of  the  generating  aircraft.  The  vortex  strength, 
quantified  by  the  circulation,  T,  is 


L  _  2 VJbC^ 

PKA  nAR 


(3.5. 1.1) 
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The  vortex  flow  velocity  on  aircraft  flying  in  formation  is  dependent  on 
circulation,  the  distance  from  the  vortex  to  the  trailing  aircraft  wing,  r,  and  the  radius  of 
the  vortex  core,  rc.  The  vortex  core  radius  is  defined  as  the  distance  from  the  vortex 
center  to  the  point  of  maximum  velocity  and  is  a  function  of  the  kinematic  viscosity  of 
the  fluid.  The  Burnham  profile  for  velocity  is, 


V  =  — 

vor'  In 


\ 


r 2  +  r 2 
V  c  J 


(3.5.1.2) 


Applying  the  superposition  principle  to  the  trailing  vortices  and  neglecting  the 
filament  attached  to  the  wing  of  L,  the  vortex  flow  velocity  is  the  difference  between  the 
flow  from  the  near  vortex  and  the  flow  from  the  far  vortex.  This  relationship  is 
illustrated  in  Figure  5.  However,  the  complication  arises  in  the  calculation  of  ri  and 
the  distance  from  the  vortex  filaments  to  a  point  on  the  wing,  because  now  the  wing 
separation  and  attitude  must  be  expressed  in  the  L  reference  frame.  By  using  vector 
addition  and  axes  transformations,  all  distances  are  resolved  into  the  L  rotating  reference 
frame.  Distances  from  the  vortices  to  a  point  on  the  trailing  aircraft  can  be  found  from 
geometry. 
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(3.5. 1.3) 


/ 

0 

X 

'o' 

r  y\  ^ 

XL 

-Kl 

-CLCr 

y 

~cLc' 

yL  • 

V 

0 

z 

0_ 

/ 

A. 

3-13 


The  angle  associated  with  each  distance,  is  the  inverse  cosine  of  the  y-direction 
component  of  the  corresponding  distance  divided  by  the  whole  distance.  Note  that  £  is 


the  integration  variable  used  to  represent  the  distance  along  the  wing  span. 


0l  =  cos 


-i 


v  1  J 


(3.5.1.4) 


62  =cos  11  '2 


V  2  ) 


The  velocities  caused  by  the  vortices  can  be  superposed  to  give  a  total  average 
velocity  vector  in  the  L  rotating  reference  frame, 


u 
v 
w 

L^J  L 


0 

Korn  sin#!  ~Vvorl2sin62 
Korn  COS01  -Vmr,2COsd2 


<yL> 

Zl 


Substituting  for  the  vortex  velocities,  Vvort,  with  the  Burnham  profile  expression  found  in 
equation  (3. 5. 1.2)  and  the  angle  relationships  from  equation  (3. 5. 1.4)  creates, 
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This  can  be  further  modified  by  substituting  the  relationship  for  circulation,  canceling 
appropriate  terms,  and  inserting  the  applicable  distance  relationships.  The  xL-component 
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of  distance  is  set  to  zero  because  of  the  assumption  of  constant  vortex  strength  of  the 
vortex  filament  in  the  XL-direction. 
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Ll 


2  npVJvL 
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A. 

y  l 

z, 


(3.5.1.5) 


Using  the  transformation  matrices  the  vortex  related  velocity  vector  can  be  transformed 
from  the  L  rotating  reference  frame  into  W. 
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(3.5.1.6) 
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3.5.2  Lift  and  Drag  Analysis 


ZW  + 


Xw 


y w  out  of  page 


(View  A) 


Figure  6.  Lift  Vector  Rotation  on  W  Wing. 


View  A  in  Figure  6  shows  the  wing  in  ordinary  flight,  meaning  a  single  aircraft 
with  no  vortex  interaction.  View  B  shows  the  rotation  of  the  lift  vector  due  to  the  action 
of  the  L  vortex  filament  and  the  new  lift,  L',  in  the  direction  of  the  original  lift  vector.  L' 
and  D'  are  in  the  direction  of  the  original  lift  and  drag,  but  increased  and  decreased  in 
magnitude,  respectively.  The  reduction  of  induced  drag  is  the  essence  of  close  formation 
flight.  The  change  in  the  angle  of  attack  as  shown  in  Figure  6  can  be  found  using 
trigonometry.  The  assumption  of  a  small  incidence  angle,  E,  is  substantiated  in  [15], 
where  the  research  found  £  ranged  from  two  to  four  degrees.  Applying  small  angle 
assumptions,  the  inverse  tangent  of  a  ratio  is  approximately  equal  to  the  ratio. 
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£=tan  l{w!Vm)  ~  w/V0 


The  next  step  is  to  derive  expressions  for  the  change  in  lift  and  drag  from  the 
original  values  to  the  new  values.  This  is  done  by  starting  with  the  definition  of  lift  as 
being  the  product  of  the  lift  curve  slope,  La,  and  the  angle  of  attack,  a.  The  subscript  o 
refers  to  the  original  lift  or  drag. 


K  =  Ka 


r 

L  = 

f  D ^ 

o 

0 

J; 

Similarly  the  rotated  lift  and  drag,  subscript  a,  can  be  written  as, 

La=La(cc  +  e) 


(3.5.2. 1) 


£>  = 


La  = 

<LJ 

a 

lLJ 

La(a  +  s) 


The  small  angle  approximation  is  also  used  to  simplify  the  trigonometric 
expressions  for  the  new  lift  and  drag  associated  with  the  flight  path  direction. 


L'  =  La  cos  £  +  D a  sin  £  ~  La  +  Da  £ 


D'=  Da  cos  £  -  La  sin  £  «  Da  -  Lae 


(3. 5.2. 2) 


Substituting  the  expressions  from  equation  (3.5.2. 1)  into  the  L'  equation  shown  above 
yields, 


f  D ' 


L'=La(cc  +  e)+  —  La(a  +  e)s 


\  L  j 
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Multiplying  and  collecting  like  terms  allows  the  new  lift,  L'  to  be  written  as, 


L'=Lacc  +  L0 


e  + 


ne  + 


f  ZY 


\ 


For  the  small  values  attributable  to  e,  it  is  a  valid  approximation  to  say  that  the  original 
drag  to  lift  ratio  is  equal  to  the  rotated  drag  to  lift  ratio,  and  that  the  lift  curve  slope 
remains  constant.  The  new  lift,  L',  can  be  written  in  terms  of  the  original  lift,  Lq,  by 
making  substitutions  and  neglecting  the  very  small  terms. 


L'=L0+Lae 

Therefore  the  change  in  lift  resulting  from  vortex  interaction  can  be  written  as  a  function 
of  lift  curve  slope,  freestream  velocity,  and  the  upwash  velocity,  w. 


A L  =  La 


(  w  \ 


(3. 5.2.3) 


A  derivation  similar  to  the  one  for  the  change  in  lift  can  be  conducted  to  find  the 
new  drag,  D'.  Again  the  difference  between  the  original  drag  and  the  new  drag  can  be 
found  starting  with  substituting  the  relationships  for  Da  and  La  from  equation  (3.5.2. 1) 
into  the  expression  for  D'  from  (3. 5. 2.2). 


D'= 


U, 


La(a  +  e)-L«(«  +  e)e 


Again  multiplying  and  collecting  like  terms  allows  D'  to  be  written  as, 


f  dY 

( 

9 

f  D) 

\ 

Laa-Ln 

GCE  +  E  - 

£ 

UJ 

cc  a 

V 

(L) 

/ 
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As  in  the  new  lift  derivation,  assumptions  are  made  about  constant  drag  to  lift  ratios  and 
small  angles  which  yield, 

D'=  D0-Laae 


As  with  the  change  in  lift,  the  change  in  drag  caused  by  close  formation  flight  can  be 
written  as  a  function  of  lift  curve  slope,  angle  of  attack,  freestream  velocity,  and  the 
upwash  velocity,  w. 


AD  =  -Laa 


w 

V~ 


(3. 5.2.4) 


It  is  vital  to  recognize  that  the  z-component  of  velocity,  w,  used  in  this  derivation 
is  in  the  W  rotating  reference  frame. 


3.5.3  Upwash  Velocity  Equations 

From  equations  (3. 5.2.3)  and  (3.5.2.4),  the  only  information  remaining  to  be 
found  in  the  calculation  of  the  change  in  lift  and  drag  is  the  of  upwash,  w,  across  the  W 
wing.  This  is  accomplished  by  first  finding  the  average  velocity  vector  components  in 
the  L  frame  of  reference  and  then  transforming  the  vector  into  the  W  frame.  The 
integration  is  conducted  as  though  the  wing  were  a  lifting  line  equal  in  distance  to  the 
span  between  the  trailing  vortices,  bv,  and  is  non-dimensionalized  by  dividing  the 
integrand  by  bv. 
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■w'2 

f  La 

f  \ 

w  1 

dy >1 

J  a 

VlV  t  2 

l^j 

bvW  /  2 

(A D)w  —  J  ”  Lacc 

Kw  1 2 


rwYdy' 


yvW 


By  substituting  the  relationship  for  w  in  the  integrand,  which  is  the  vortex  velocity 
component  in  the  zw-direction,  with  the  velocity  expression  from  equation  (3.5. 1.6)  and 
removing  the  constants  from  the  integrand  yields  the  following  expressions: 


(AL)W  = 


Vb 


vW 


Ll 


2npV„bvL  j 


(A D)w=~L°a 


Vb 


vW 


2  npVJb 


vL  J 


/~1  W  /  _ 

C,  CLa 


W  I  _ 

C,  CLa 


Where  the  function  a  is  defined  as. 


a  = 


W  f  2 

f 

(rM) 

f 

J 

\w 1 2 

kA 

W^j 

w/ 2 

f 

-ifo}  \ 

J 

W/2 

r* 

k 


y  l 

lA  J 


with  the  expressions  for  ri  and  r2  found  in  equation  (3. 5. 1.3).  Integrating  and  collecting 
terms  yields, 


cr(xj  =  0 


(3.5.3. 1) 
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^(>0  = 


-z 

n  tan 

bvC-bv+2Cy  ' 

■yjc2x2  +  C2z2  +  r2 

U  UU 1 

2 Jc2x2  +C2z2  +r?  j 

a  tan 


a  tan 


v2^/cV+CV+J  ; 

-bvC  +  bv  +  2Cy 
2  ,Jc2x2+C2z2+r? 


-  a  tan 


(  _^b£_-K+2Cy_^ 
2  JcV  +  CV+rc! 


+ 


(3. 5. 3. 2) 


^(zJ 


_1_ 

2C 


In 


r  Z>v(l  +  2C  +  C2)+  4C6„;y(l  +  C)  +  4C2(x2  +  y2  +  z2)+4/;2  " 
^(l  -  2C  +  C2)+4Cfcl,y(-l  +  c)+  4C2(x2  +  y2  +  z2)+4rc2  j 


* 


^(l  +  2C  +  C2)+4Cfevy(-l-C)  +  4C2(x2  +  y2  +z2)+4rc2^ 
fcv(l  -  2C  +  C2)+4CZ?vy(-l  +  C)  +  4C2(y2  +  y2  +  z2)+  4rc2  ^ 


(3.5. 3. 3) 


where  C  is  short  hand  for  the  coordinate  transformation, 

C  =  C’fCl 

Note  that  when  the  two  rotating  reference  frames  are  aligned,  these  three 
equations  are  mathematically  equal  to  the  ones  found  in  both  [8]  and  [10]. 

3.5.4  Vortex  Effects  on  the  Vertical  Tail 

Recalling  the  requirement  to  fly  with  no  sideslip,  the  only  side  force,  Yw,  on  the 
vertical  tail  of  the  wing  aircraft  is  caused  by  the  component  of  vortex  flow  in  the  y- 
direction  of  the  W  rotating  reference  frame.  The  analysis  for  the  side  force  follows  the 
same  steps  as  wing  analysis,  the  primary  difference  coming  in  the  distance  vector 
calculation. 
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Figure  7.  Vertical  Tail  of  W  in  Vortex  Interaction. 


As  shown  in  Figure  7,  the  distance  from  the  vortices  to  a  location  on  the  vertical 
tail  can  be  written  as  vectors  in  the  associated  rotating  reference  frame,  which  can  then  be 
transformed  so  all  distances  are  in  the  L  frame.  The  only  difference  between  the  wing 
and  vertical  tail  separation  is  the  integration  variable,  %  that  now  is  in  the  z  direction  of 
the  W  reference  frame. 


(3.5.4. 1) 
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The  force  and  drag  analysis  on  the  vertical  tail  follows  a  similar  pattern  as  for  the 
wing  with  all  the  same  assumptions  and  conditions,  except  that  the  initial  angle  of  attack 
is  zero.  The  change  in  side  force  is  the  parameter  of  interest  and  has  the  same  result  as 
the  change  in  lift  on  the  wing.  Figure  8  shows  the  velocity  vectors  and  forces  on  the 
vertical  tail. 


v 


Figure  8.  Interaction  Forces  on  the  Vertical  Tail  (viewed  from  above). 

The  change  in  side  force  on  the  vertical  tail,  AY,  can  be  written, 

A  Y  =  Lajvl 

Again,  the  form  of  the  vortex  interaction  equation  remains  the  same  as  for  the 
wing,  with  the  function  avt  and  lift  curve  slope  being  subscripted  to  reflect  the  vertical  tail 
values,  as  opposed  to  those  of  the  wing. 
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K  (  h 


v.K 


2npVJvL 


s~iW  /~t I  — . 

c,  CLov 


Where  orvt  is  defined  as. 


^v,  = 


0 

V  ,  O  f 

f  riw  j£z,  J 
o[hw2+^-2  , 

*'y  r„{yJl 


Mvf _ 

L0Vriw2+rc2  J 


2  2 
/2w  +r, 

r  2  +  r2 

A2v/  t  'c 


XL 

dE, 

A 

y  l 

{ 

ni> 

r- 

k 

The  limits  of  the  integration  are  from  zero,  the  longitudinal  centerline  of  W,  to  bvt, 
the  height  of  the  vertical  tail.  The  results  of  the  integration  are, 

ct(xJ  =  0  (3. 5.4. 2) 


2y 


^4C2x2+{2Cy-bvJ+4r2 


( 


a  tan 


2Cz 


a  tan 


4C2x2  +  (2  Cy  +  bvL  )2  +  4rc2 


a  tan 


( 


a  tan 


^4C1x2 

+  (2Cy-b,J 

+  4r2 

2C(b„  +  z) 

^4  CV 

+  (20*  -  bvL  )2 

+  4r2 

2Cz 

^4  C2x2 

+  (2Cy  +  b,J 

+  4r2 

2  C(-b„+z) 

J* 4C2x 2 

+  {2Cy  +  bvL  f 

+  4  r2 

j  j 


+ 


yj 


(3.54.3) 
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4  Formation  Flight  Controller 


4.1  Controller  Design 

The  foundation  for  the  controller  design  is  the  PID  control  action  on  each  of  the 
three  control  channels.  For  the  primary  control  law  thrust  control  is  affected  in  the  x 
channel,  which  is  to  say  that  any  error  in  the  commanded  and  the  actual  value  of  x 
separation  distance,  results  in  a  change  in  thrust.  In  a  similar  manner,  lift  is  the  control 
variable  in  the  z  channel  and  roll  rate  is  the  control  variable  in  the  y  channel.  Eqs.  (4.1.1) 
-  (4.1.3)  show  the  actual  control  law  used  for  each  of  the  three  control  variables:  Tw,  Lw, 
and  pw,  respectively.  Note  that  the  dynamics  associated  with  lateral  maneuvers  required 
the  inclusion  of  a  second  derivative  on  separation  error  and  proportional  feedback  of  the 
heading  error  for  the  roll  rate  control  variable. 


Tw  —  Kj,  e %  +  K j,  ex  +  K j,  \e ^dt 


L\v  =  Kipez  +*lK  +Kl j  ^ezdt 


pw  =  K  e  +  K  ew  +  K  e  +  K  e  +  K  j  e  dt 
Fw  Pp  y  Wep  Vt  Pd  y  Pdd  y  Pi  y 


(4.1.1) 

(4.1.2) 

(4.1.3) 


Where  the  error  signals  are  the  difference  between  the  commanded  and  actual  formation 
geometry  separation  distances. 


(4.1.4) 
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The  integrals  of  the  error  signals  are  created  by  augmenting  the  formation 
equations  of  motion.  The  derivative  control  elements  are  created  by  taking  the  difference 
in  the  separation  distances  for  a  single  time  step  and  dividing  by  the  magnitude  of  the 
time  step.  As  noted  above,  control  in  the  lateral  channel  required  additional  control 
elements  in  the  control  law. 

In  an  attempt  to  more  closely  approximate  actual  pilot  control  inputs,  an 
alternative  control  law  was  devised.  For  this  new  control  law,  the  need  for  W  to  climb, 
according  to  a  separation  error  in  the  x  direction,  results  in  the  controller  responding  with 
increased  thrust,  whereas  the  need  for  W  to  increase  flight  speed  causes  the  lift  to 
decrease.  In  other  words,  thrust  is  now  the  control  variable  in  the  z  channel  and  lift  is  the 
control  variable  in  the  x  channel.  Roll  rate  continues  as  the  control  variable  in  the  y 
channel.  Another  attempt  to  simulate  the  "seat  of  the  pants"  feedback  used  by  pilots  is 
the  inclusion  of  acceleration  feedback  in  the  alternative  control  law.  Equations  for  the 
new  control  law  are  shown  below. 

Tw=KTe  +KTe  +K  e  +K  je  dt  (4.1.5) 

I  p  Z  l  D  Z  *  DD  Z  1 1  X 

LW=K  e+K  e+K  e  +  K ,  je  dt  (4.1.6) 

A  Up  A  UDD  A  U{  4 

pw=K  e  +K  em  +K  e+K  e+K  jedt  (4.1.7) 

yw  pP  y  weP  v '  pD  y  Pdd  y  Pi  y 
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In  the  expressions  for  the  control  law  above,  pw  is  called  the  roll  rate  control 
variable.  The  actual  roll  rate,  Pw,  can  be  modified  from  what  is  shown  in  equation 
(3.2.6),  by  virtue  of  a  dynamic  inversion  argument,  to  be  written  as. 


Pw  =  K  e  +  Kt 
W  Pp  y  Vep 


ew  +K  ev+K  ev+K  fe  dt  + 
Pd  y  Pdd  y  Pi  y 


tan  Yw  {Lw  sin  (j)w  -  Yw  cos  <j)w ) 


mw  Vw 


4.2  Formation  Flight  Turns 

The  formation  flight  control  system  is  such  that  right  and  left  turns  have 
dissimilar  dynamics,  dependent  upon  W's  position.  This  is  similar  to  running  on  a  track, 
where  the  runner  on  the  inside  lane  has  to  turn  sharper  (on  a  shorter  path)  to  stay  in  his 
lane.  Moreover,  pilots  fly  three  different  formations  during  turns.  In  the  first  type  of 
formation  turn,  L  remains  in  W's  x-y  plane.  This  is  the  method  used  for  the  design  of  the 
automatic  formation  flight  control  system.  In  a  "wingtip"  formation,  W  remains  in  L's  x- 
y  plane  throughout  the  turn.  In  this  case,  the  error  signals  in  Eq.  (4.1.4)  are  transformed 
into  the  inertial  frame  of  reference  according  to 


1 

1 _ 
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1 

X 
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II 

yc-y 

_ez_ 

- ! 

N? 

1 

_ 1 

prior  to  the  application  of  the  PED  control  law  (Eqs.  (4.1.1)  -  (4.1.3)).  In  the  third  mode, 
called  a  "route"  formation,  W  and  L  remain  at  the  same  inertial  x-y  plane  during  the  turn. 
This  type  of  formation  would  be  accomplished  in  the  simulation  by  transforming  the 
error  signals  in  Eq.  (4.1.4)  into  the  L  reference  frame  according  to 
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again  prior  to  the  application  of  the  PID  control  law. 
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5  Formation  Flight  Simulation  Results 


5.1  Simulation  Aircraft 

The  formation  flight  control  system  response  to  maneuvers  flown  by  the  leader 
and  to  commanded  changes  in  formation  geometry  (xc,  yc,  and  zc)  was  simulated  using 
the  equations  of  motion  for  the  two  aircraft  system  shown  in  equations  (3.4.2)  -  (3.4.8). 
Both  aircraft  used  in  the  formation  flight  simulations  were  representative  of  F-16  class 
aircraft.  In  the  simulations,  the  formation  is  flying  at  an  altitude  of  15,000  meters  and  at 
an  airspeed  of  0.85  Mach,  or  approximately  251.5  meters  per  second.  The  simulations 
were  conducted  with  the  assumption  of  incompressible  flow.  Assumptions  were  made  as 
to  the  values  of  the  constants  in  the  drag  polar  equation  (shown  in  Table  1)  in  order  to 
determine  the  thrust  required  for  steady  level  flight.  Naturally,  the  lift  balanced  the 
weight  during  trimmed  flight.  The  following  table  shows  the  aircraft  characteristics, 
obtained  from  [8],  used  in  the  simulations. 


Table  1  F-16  Class  Aircraft  Specifications 


Wing  Area 

S 

27.87  nr 

Wing  Span 

b 

9.14  m 

Aspect  Ratio 

AR 

3 

Lift  Curve  Slope 

La 

5.3  per  rad 

Coefficient  of  Zero-Lift  Drag 

CDo 

0.015  1 

Drag  Polar  Constant 

K 

0.02  ! 

Tail  Area 

Svt 

5.086  n? 

Tail  Height 

bvt 

3.05  m 

Tail  Lift  Curve  Slope 

K 

5.3  per  rad 

Efficiency  Factor 

y\ 

0.95 

Mass 

m 

11336.4  m 
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Given  the  wing  span  annotated  in  Table  1  and  referring  to  the  formation  geometry 
discussion  in  Section  2.2.3,  the  formation  geometry  used  in  the  simulations,  measured  in 
the  W  reference  frame,  was  xc  =  27  meters,  zc  =  0  meters,  and  yc  =  7  meters. 


5.2  Displacement  Recovery  Maneuvers 

Initially  the  formation  flight  control  system  was  exercised  by  commanding  W  to 
perform  a  repositioning  after  being  displaced  by  one  meter  in  each  of  the  three  Cartesian 
axes  directions.  The  control  system  gains  were  selected  to  ensure  the  aircraft  would 
properly  return  to  the  commanded  position.  The  simulation  results  of  the  control  system 
performance  are  shown  by  way  of  time  histories  for  each  maneuver.  For  several 
maneuvers,  the  alternative  control  system  response  is  also  shown  for  comparative 
purposes. 
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meters 


Figure  9  displays  the  control  system  response  to  a  one  meter  step  input  (x  =  xc+l, 


y  -  yc+l,  and  z  =  zc+l)  to  the  separation  distances. 
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Figure  9. 


W  displaced  from  the  optimal  position. 


In  Figure  10,  the  alternative  control  law  was  used  to  reposition  W  under  the  same 
initial  conditions  as  in  the  previous  figure.  Comparing  the  two  system  responses,  the 
alternative  control  law  takes  longer  to  recover  and  has  larger  control  inputs. 


X  Separation  x  1 0s  W  Lift 


0  20  40  60  0  20  40  60 


Z  Separation  W  Thrust 


Figure  10.  W  displaced  from  optimal  position:  alternate  control  law 
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In  Figure  11,  the  one  meter  step  input  has  been  changed  (x  =  xc-l,  y  =  yc-l,  and 
z  =  zc-l)  to  ensure  the  control  system  is  versatile  enough  to  reposition  from  a  different 


initial  displacement.  Again  the  controller  restored  W  to  the  correct  formation  position. 
X  Separation 


Y  Separation 


Thrust 


10 


20 


30 


Z  Separation  x  -jq5  Lift 


Figure  11.  W  displaced  away  from  the  optimal  position. 
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Because  of  symmetry,  the  formation  can  be  aligned  on  either  the  right  or  left  wing 
of  the  leader.  In  the  formation  flight  control  system,  this  is  accomplished  by  simply 
changing  the  sign  on  the  formation  geometry  value  of  commanded  y  separation,  yc.  For 
this  simulation  positive  seven  became  negative  seven  and  W  was  displaced  from  this  new 
commanded  formation  position.  The  controller  response  shown  in  Figure  12  is  evidence 
that  the  formation  can  be  located  on  either  side  of  L. 


X  Separation 


Thrust 


Y  Separation 


Bank  Angle 


Z  Separation  x  105  Lift 


Figure  12.  W  displaced  from  the  new  optimal  position. 
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5.3  Velocity  Change  Maneuvers 


The  next  step  was  to  fly  L  through  several  maneuvers  and  adjust  the  gains  to  the 
levels  required  to  maintain  the  formation.  Three  different  maneuvers  were  explored:  L 
accelerating  to  a  new  velocity,  L  changing  altitude,  and  L  changing  heading.  It  is 
important  to  note  that  the  requirement  of  no  sideslip  mandates  that  all  heading  changes  be 
accomplished  by  rolling  the  aircraft.  While  rolling  the  aircraft,  altitude  and  velocity  are 
maintained  by  adjusting  lift  and  thrust.  In  other  words,  making  the  heading  change  is  a 
kinematically  coupled  maneuver,  where  adjustments  to  any  control  channel  affects  the 
other  channels.  Many  trials  were  made  to  find  suitable  gains  to  accomplish  the  heading 
change  maneuver.  These  gains  were  then  reapplied  to  the  other  maneuvers  to  ensure  the 
system  would  still  be  able  to  recover  the  formation  geometry.  The  first  simulations,  see 
Figure  9  -  Figure  12,  involved  L  flying  straight  and  level  and  W  recovering  from  an 
initial  displacement.  The  next  series  of  simulation  results  have  the  leader  fly  various 
maneuvers  that  could  be  expected  during  a  normal  flight.  The  simplest  of  these 
maneuvers  was  a  change  in  velocity,  both  positive  and  negative. 

Figure  13  shows  how  W  responded  when  L's  velocity  was  increased  by  10  meters 
per  second.  The  velocity  lag  in  W  is  minimal,  and  the  thrust  increases  in  conjunction 
with  the  x  separation  distance  and  settles  at  a  new  steady  state  value  which  is  slightly 
greater  than  the  original  trim  point,  as  expected. 
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Figure  13.  W's  response  to  L  velocity  increase  maneuver. 


For  comparison  purposes,  the  velocity  increase  maneuver  is  repeated  using  the 
alternate  control  law.  Again,  the  indirect  nature  of  the  alternative  control  law  causes  the 
second  controller  to  respond  differently  to  system  disturbances.  In  particular,  note  in 
Figure  14  the  difference  in  the  x  separation  time  histories  and  the  significant  amount  of 
lift  adjustment  made  by  the  second  controller  when  changing  velocity,  whereas  for  the 
primary  control  system,  the  change  in  lift  during  the  velocity  change  maneuver  was 
almost  negligible. 


5-8 


L  Velocity 


W  Velocity 


X  Separation 


Figure  14.  W's  response  to  L  velocity  increase  maneuver:  alternate  control  law. 


In  the  next  velocity  change  maneuver,  shown  in  Figure  15,  L  decreased  airspeed 
and  the  automatic  controller  maintained  the  formation.  Note  that  the  value  of  peak  thrust 
is  actually  negative.  If  the  model  were  of  a  large  transport  airplane,  this  could  be  thought 
of  as  engine  thrust  reversing,  but  for  an  F-16,  negative  thrust  would  come  from  the 
increased  drag  of  air  brakes  being  deployed  as  the  throttles  were  pulled  back. 


5-9 


L  Velocity 


W  Velocity 


X  Separation  W  Thrust 


Figure  15.  W's  response  to  L  velocity  decrease  maneuver. 


5.4  Altitude  Change  Maneuvers 

Another  pitch  plane  maneuver  is  the  formation  altitude  increase.  In  the  figure 
below,  the  pitch  angle  of  W  changes  from  zero  degrees  to  negative  five  degrees  and  back 
to  zero  as  W  climbs  to  the  new  altitude.  This  is  due  to  the  definition  of  positive  pitch 
angle  in  the  wind  axes  system,  as  illustrated  in  Figure  3.  L  is  flown  at  a  constant 
velocity,  which  causes  the  control  system  to  increase  W's  thrust  to  maintain  proper  x 
separation  throughout  the  climbing  portion  of  the  maneuver. 
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degrees  meters  meters 


Note  in  Figure  16  that  the  z  separation  shows  the  control  system  lag  as  W  adjusts  to  the 
leader  climbing  and  then  leveling  at  the  new  altitude. 


X  Separation  W  Thrust 


W  Pitch  Angle  x  ^q4  Formation  Altitude  Change 


seconds  seconds 


Figure  16.  W's  response  to  L  altitude  increase  maneuver. 
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To  compare  the  second  controller  to  the  one  of  primary  interest,  the  altitude 
increase  maneuver  was  repeated.  Note  in  the  figure  below  that  the  peak  values  of  thrust 
and  lift  are  nearly  identical  to  those  of  the  primary  control  law.  However,  the  transients 
in  the  separation  distances  are  dramatically  larger  using  the  second  controller. 


X  Separation  x  105 


Z  Separation 


W  Thrust 


W  Pitch  Angle  x  104  Altitude  Change 


Figure  17.  W's  response  to  L  altitude  increase  maneuver:  alternate  control  law. 
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degrees  meters  meters 


The  controller  time  histories  shown  in  Figure  18  display  the  control  system  taking 
W  to  a  lower  altitude  in  response  to  L  descending.  Note  that  for  this  maneuver,  thrust 
becomes  negative,  similar  to  the  velocity  decrease  maneuver. 


X  Separation  W  Thrust 


Z  Separation  x  105  W  Lift 


W  Pitch  Angle 


xIO4  Formation  Altitude  Change 


Figure  18.  W's  response  to  L  altitude  decrease  maneuver. 
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5.5  Heading  Change  Maneuvers 


The  kinematic  coupling  that  occurs  during  turns  causes  the  heading  change 
maneuvers  to  be  more  complex  for  the  system  to  control.  The  first  turning  maneuver  L 
conducted  was  a  positive  heading  change  while  maintaining  constant  altitude  and 
velocity.  Due  to  the  definition  of  the  wind  axes  system  (Figure  3),  a  positive  heading 
change  is  produced  by  turning  to  the  left.  In  this  maneuver,  the  formation  geometry  has 
W  positioned  to  the  right  of  L,  and  therefore  turning  along  a  wider  track.  Figure  19 
shows  roll  angle  and  lift  of  L  as  it  makes  the  heading  change. 


Figure  19.  L  positive  heading  change  maneuver. 


The  next  series  of  figures  show  the  results  of  W's  control  system  responses  on  its 
position  and  attitude  relative  to  L.  The  separation  distances  measured  in  the  inertial 
reference  frame  have  been  included  to  show  the  differences  that  occurred  during  this 
maneuver.  For  each  control  channel,  the  alternate  controller  response  to  the  same 
disturbance  signal  is  shown. 
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The  velocity  increase  in  Figure  20  shows  the  effect  of  the  larger  turn  radius  flown 
by  W.  In  order  to  maintain  the  formation  W  has  to  travel  at  a  greater  velocity.  Note  the 
difference  in  the  final  value  of  the  inertial  separation  distance  compared  to  the  W 
reference  frame  separation  distance. 


X  Separation 


X  Separation  (Inertial) 


W  Thrust 


W  Velocity 


Figure  20.  W's  response  to  positive  L  heading  change:  x  axis. 


Comparing  the  x  channel  response  of  the  primary  control  system,  shown  above, 
with  the  alternate,  shown  in  Figure  21,  the  x  separation  actually  decreases  during  the  turn, 
which  is  validated  by  the  corresponding  increase  in  velocity.  The  lift  commanded  by  the 
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primary  controller,  see  Figure  24,  is  only  slightly  less  than  that  commanded  by  the 


alternate  controller. 

X  Separation  x  105  W  Lift 


X  Separation  (Inertial)  W  Velocity 


Figure  21.  W's  response  to  positive  L  heading  change:  alternate  control  law:  x  axis. 


The  y  separation  distances  shown  in  Figure  22  are  the  results  caused  by  L  banking 
to  the  left  and  flying  away  from  W  initially  and  then  rolling  back  to  level  flight  and 
having  W  overshoot  the  commanded  y  separation  distance.  Comparing  the  roll  angle  of 
W  (Figure  22)  with  that  of  L  (Figure  19),  the  time  histories  are  nearly  identical,  with  only 
a  slight  overshoot  by  the  control  system,  as  expected.  Again  the  two  separation  distances 
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are  dramatically  different.  These  distances  are  the  Cartesian  components  of  the  total 
separation  distance.  The  magnitude  of  the  formation  separation  distance  is  equal  in  all 
reference  frames. 


Y  Separation  W  Roll  Angle 


Y  Separation  (Inertial)  W  Heading  Angle 


Figure  22.  W's  response  to  positive  L  heading  change:  y  axis. 


The  alternate  control  law  has  very  little  effect  on  the  y  channel  response,  except 
that  the  system  is  slower  and  the  separation  distance  transients,  shown  in  Figure  23,  are 
greater  than  for  the  primary  control  system.  Although  the  response  is  different,  both 
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separation  distances  have  the  same  steady  state  value  as  for  the  primary  control  system 


response. 

Y  Separation  W  Roll  Angle 


Y  Separation  (Inertial)  W  Heading  Angle 


Figure  23.  W's  response  to  positive  L  heading  change:  alternate  control  law:  y  axis. 

The  z  channel  controller  response  shown  in  Figure  24  provided  a  dramatic 
increase  in  lift,  as  would  be  expected  to  maintain  altitude  during  a  coordinated  turn.  Note 
the  large  disagreement  between  the  two  reference  frame  measurements  of  the  z  separation 
distance  transients.  The  inertial  frame  shows  a  much  larger  value  of  z  separation, 
indicating  the  strong  effect  of  roll  angle  on  measured  distances.  The  comparatively  small 
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value  of  separation  distance  in  the  W  frame  is  the  result  of  the  formation  turn  mode  used 
in  the  control  system,  meaning  that  the  controller  is  trying  to  keep  L  in  W's  x-y  plane,  as 
explained  in  Section  4.2.  The  separation  distance  errors  used  in  the  controller  are 
measured  in  W's  reference  frame. 


Z  Separation 


Z  Separation  (Inertial) 


xIO5  W  Lift 


0  20  40  60 

seconds 


Figure  24.  W's  response  to  positive  L  heading  change:  z  axis. 


To  compare  commanded  thrust  of  the  two  control  laws  refer  to  Figure  20  for  the 
primary  and  Figure  25  for  the  alternate.  The  alternate  control  system  requires  larger 
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thrust  variations  to  maintain  the  formation.  As  seen  in  the  x  and  y  axes,  the  separation 
distance  transients  are  larger  for  the  alternate  controller  than  for  the  primary. 


Z  Separation 


W  Thrust 


Z  Separation  (Inertial)  W  Pitch  Angle 


Figure  25.  W's  response  to  positive  L  heading  change:  alternate  control  law:  z  axis. 


Section  4.2  discusses  optional  geometry  as  the  formation  conducts  turning 
maneuvers.  The  principal  approach  in  the  simulations  was  to  have  the  controller 
maintain  L  in  W's  x-y  plane,  which  is  to  say,  a  z  separation  distance  equal  to  zero.  Figure 
26  shows  a  "route"  formation  geometry,  where  the  controller  maintains  L  and  W  flying  in 
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meters 


the  same  inertial  x-y  plane.  Notice  that  the  separation  distance  in  the  inertial  reference 
frame  is  very  small  compared  to  the  W  reference  frame. 

Z  Separation 


Z  Separation  (Inertial)  W  Pitch  Angle 


seconds  seconds 


Figure  26.  W's  response  to  positive  L  heading  change:  route  formation:  z  axis. 
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The  next  maneuver  run  in  the  simulation  was  a  heading  change  to  the  right  rather 
than  left,  or  in  other  words,  L  is  turning  to  a  negative  heading  angle,  \|/,  as  shown  in  the 
figure  below. 


L  Roll  Angle  x  105  L  Lift  L  Heading  Angle 


Figure  27.  L  negative  heading  change  maneuver. 
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When  L  turns  to  the  right,  W  is  now  flying  on  the  short  track.  Therefore,  the 
controller  must  reduce  W's  velocity  to  prevent  overshooting  the  leader,  as  shown  in 
Figure  28.  As  noted  previously,  the  x  component  of  the  inertial  separation  distance 
reflects  the  geometry  of  the  formation  from  a  different  perspective,  but  does  not  change 
the  magnitude  of  the  separation  distance. 


X  Separation  W  Thrust 


X  Separation  (Inertial) 


W  Velocity 


20  40 

seconds 


60 


Figure  28.  W's  response  to  negative  L  heading  change:  x  axis. 
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It  is  evident  from  Figure  29  that  the  control  system  is  capable  of  compensating  for 
the  leader  turning  left  or  right.  Again,  the  roll  angle  of  W  matches  L  almost  perfectly, 
with  only  a  slight  overshoot.  The  steady  state  value  of  inertial  separation  distance 
decreased,  as  expected. 


Y  Separation  (Inertial) 


W  Heading  Angle 


seconds 


seconds 


Figure  29.  W's  response  to  negative  L  heading  change:  y  axis. 
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meters 


As  for  the  positive  heading  change,  Figure  30  shows  a  dramatic  difference  in  the  z 
separation  of  the  rotating  and  inertial  reference  frame,  highlighting  the  dependency  of 
separation  distances  on  reference  frame  orientation. 


Z  Separation  x  105  W  Lift 


seconds  seconds 

Figure  30.  W's  response  to  negative  L  heading  change:  z  axis. 


5-25 


meters 


5.6  Climbing  Turn  Maneuver 

Most  complex  maneuver  conducted  in  the  simulation  was  the  climbing  turn. 
Figure  31  shows  the  three  dimensional  flight  path  of  the  formation  as  it  climbed  to  a  new 
altitude  and  turned  left. 


x  104 


Flight  Path  From  Behind 


Figure  31.  Flight  path  in  climbing  turn  maneuver:  inertial  reference  frame. 
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To  produce  the  altitude  and  heading  change  in  the  formation,  the  disturbance 
signals  shown  below  were  created  by  L.  The  amount  of  heading  change  and  the  distance 
climbed  for  the  combined  maneuver  is  less  than  the  individual  maneuvers,  because  the 
more  aggressive  maneuvers,  when  combined  together,  caused  the  controller  to  be 
unstable. 


L  Roll  Angle 


xIO5  LLift 


L  Pitch  Angle 


0  20  40  60 

seconds 


Figure  32.  L  simultaneous  altitude  and  heading  change  maneuver. 


The  next  three  figures  show  the  formation  flight  control  system  performance  and 
the  W  response,  in  the  three  Cartesian  axes. 

Comparing  the  thrust  required  to  perform  the  straight  ahead  altitude  increase 
(Figure  16)  versus  the  simultaneous  altitude  and  heading  change,  shown  below,  it  can  be 


5-27 


seen  that  turning  while  climbing  requires  significantly  more  thrust  to  maintain  the 
formation.  Because  the  heading  angle  changed  less  than  four  degrees,  the  steady  state 
value  of  inertial  separation  distances  in  both  the  x  and  y  axis  are  only  slightly  different 
than  the  W  reference  frame  distances. 


X  Separation  W  Thrust 


X  Separation  (Inertial)  VV  Velocity 


Figure  33.  W's  response  to  the  climbing  turn  maneuver:  x  axis. 
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meters  meters 


The  y  channel  of  the  control  system  was  again  able  to  match  the  L  roll  angle  to 
produce  the  correct  heading  change  with  only  slight  overshoot  in  spite  of  the  additional 
complexity  of  climbing  simultaneously,  see  Figure  34. 


Y  Separation  W  Roll  Angle 


0  20  40  60  0  20  40  60 


Y  Separation  (Inertial)  W  Heading  Angle 


seconds  seconds 

Figure  34.  W's  response  to  the  climbing  turn  maneuver:  y  axis. 
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Inasmuch  as  the  roll  angle  for  the  climbing  turn  had  a  maximum  value  of  negative 
five  degrees  (as  opposed  to  twenty  degrees  in  the  constant  altitude  turns),  Figure  35 
shows  only  a  small  amount  of  discrepancy  between  the  inertial  and  rotating  reference 
frame  values  a  z  separation  distances. 


Z  Separation  x  105  W  Lift 


seconds  seconds 

Figure  35.  W's  response  to  the  climbing  turn  maneuver:  z  axis. 
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5.7  Formation  Flight  Control  System  Gains 


Each  of  the  three  control  variables  had  proportional  control  to  be  responsive  to 
system  inputs,  integral  control  to  eliminate  steady  state  error,  and  derivative  control  to 
improve  system  stability.  The  responsiveness  was  balanced  with  the  amount  of  damping 
for  each  control  variable  individually,  and  each  control  variable  was  balanced  for  the 
overall  system.  The  final  values  of  the  control  system  gains  are  listed  in  Table  2. 


Table  2  Formation  Flight  Control  System  Gains 


Title 

Symbol 

Primary 

Alternate 

Units 

Proportional  Control  on  Thrust 

Ktp 

-1000 

-130 

N/m 

Integral  Control  on  Thrust 

Kti 

-1000 

-130 

N/(m  s) 

Derivative  Control  on  Thrust 

Kjd 

6000 

2000 

N/(m/s) 

Derivative  Control  on  Thrust 

Ktdd 

0 

5000 

N/(m/s2) 

Proportional  Control  on  Lift 

Klp 

-5000 

1300 

N/m 

Integral  Control  on  Lift 

Kli 

-7000 

1000 

N/(m  s) 

Derivative  Control  on  Lift 

Kld 

12000 

-10000 

N/(m/s) 

Derivative  Control  on  Lift 

Kldd 

0 

-1000 

N/(m/s2) 

Proportional  Control  on  Roll  Rate 

KpP 

0.008 

0.0023 

rad/(m  s) 

Proportional  Control  on  Heading  Angle 

K\|/eP 

-0.001 

0 

rad/(m  s) 

Integral  Control  on  Roll  Rate 

Kp, 

0.008 

0.0023 

rad/(m  s2) 

Derivative  Control  on  Roll  Rate 

Kpd 

-0.050 

-0.025 

rad/m 

Derivative  Control  on  Roll  Rate 

KpDD 

-0.050 

-0.025 

rad/(m/s) 

The  magnitude  of  the  gains  can  be  misleading  until  one  recognizes  the  units  associated 
each  control  element. 
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6  Conclusions  and  Recommendations 


6.1  Conclusions 

The  results  of  the  nonlinear  simulations  validate  the  concept  of  three  dimensional 
automatic  formation  flight  control.  However,  both  the  primary  and  the  alternate  control 
systems  are  very  sensitive  to  the  amplitude  of  disturbances,  which  is  to  say  that 
aggressive  L  maneuvers  can  cause  the  controllers  to  be  unstable.  Control  system 
robustness  will  become  even  more  crucial  when  the  aerodynamic  effects  of  close 
formation  flight  are  included  in  the  simulation.  Indeed  the  induced  drag  reduction 
afforded  by  close  formation  flight  requires  W  to  stay  within  +/-  5%  of  the  optimal  lateral 
separation  during  formation  maneuvers,  further  limiting  the  amplitude  of  permissible  L 
maneuvers.  Hence,  the  maneuvers  performed  by  L  must  be  relatively  benign  for  the 
controller  to  automatically  maintain  the  formation  and  the  benefit. 

6.2  Recommendations 

Further  research  on  this  project  could  be  accomplished  in  several  different  areas: 

1.  The  simulation  equations  of  motion  could  be  linearized. 

2.  The  aerodynamic  coupling  effects  of  close  formation  flight  should  be  modeled 
in  three  dimensions  and  included  in  the  simulation. 

3.  Moment  relationships  should  be  derived  and  included  in  the  simulation. 
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4.  The  current  control  law  could  be  modified  to  allow  thrust  to  be  affected  by  z 
separation  error  and  lift  to  be  affected  by  x  separation  control  error  to  reflect 
the  manner  in  which  pilots  control  their  aircraft. 

5.  The  type  of  formation  flown  during  turns  can  be  changed  to  either  of  the  other 
two  types  described  in  Section  4.1. 
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Appendix  A 


The  formation  flight  simulations  were  accomplished  using  two  Matlab  script  files, 
one  to  drive  the  simulations  and  one  that  contained  the  differentia]  equations  of  motion. 
The  Matlab  solver  routine  used  to  do  the  numerical  integration  was  called  ODE23. 


Simulation  Driver  Program 


%cf f c_sim_a . m 

global  mw  ml  g  rho  Cdow  Kw  Sw  Ll  phil  Yw  Vl  gammal; 
global  xc  yc  zc  KTI  KLI  Kpl; 

%Flight  altitude  of  15,000  m 
h=15000 ;  rho=0. 19475;  g=9.81; 

%L  aircraft  parameters 

Sl=27 . 87 ;  Llbar=111210 ;  Cdol=.015;  Kl=.02;  ml=Llbar/g;  %11336.4; 

%L  drag  relationships 
Vlbar=251 . 5 ; 
qlbar= . 5*rho*VlbarA2 ; 

Dlbar= (Cdol+Kl* (Llbar/ (qlbar*Sl) ) A2 ) *qlbar*Sl ; 

Tlbar=Dlbar ; 

%L  trim  values 

Vl=251 . 5 ;  gammal =0 ;  Tl=Tlbar;  Ll=Llbar;  phil=0;  psil=0; 

%W  aircraft  parameters 

Sw=27 . 87 ;  Lwbar=111210 ;  Cdow=.015;  Kw=.02;  mw=Lwbar/g;  %11336.4; 

%W  drag  relationships 
Vwbar=251 . 5 ;  qwbar= . 5 *rho*VwbarA2 ; 

Dwbar= (Cdow+Kw* (Lwbar/ (qwbar*Sw) ) A2) *qwbar*Sw; 

Twbar=Dwbar ; 

%W  trim  values 

uw=0;  gammaw-0 ;  phiw=0;  psiw=0;  Yw=  0 ; 

%Commanded  separation  in  the  W  reference  frame 
xc= 27;  yc=7;  zc=0; 

%X  channel  controller  gains 
KT I  =  - 1 0  0  0  ; 

KTP=-1000 ; 

KTD=6000 ; 

%Z  channel  controller  gains 
KLI=-7000 ; 

KLP=-5000 ; 

KLD=12  00  0 ; 

%Y  channel  controller  gains 
Kpl=. 0080; 

KpP= .0080; 

KpsiP=- . 0010; 

KpD=- . 050; 

KpDD=- . 050; 

%Alternate  control  law  Z  channel  controller  gains 
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%KTI=-130; 

%KTP=-130 ; 

%KTD=2000 ; 

%KTDD=5000 ; 

%  Alternate  control  law  X  channel  controller  gains 
%KLI=1000 ; 

%KLP=1300 ; 

%KLD=-10000; 

%KLDD=-1000 ; 

%  Alternate  control  law  Y  channel  controller  gains 
%KpI= . 0023 ; 

%KpP= . 0023 ; 

%KpsiP=0 ; 

%KpD=- .025; 

%KpDD=- .025; 

%Initialize  the  program  variables 

t=60 ; 
dt= . 1 ; 

n=t/dt+l;  nn=410;  nnn=0;  j=0; 
xv=0;  zv=0;  yv=0;  ya=0; 

CWI=[1  0  0 ; 0  1  0 ; 0  0  1] ; 

CLI= [ 1  0  0 ; 0  1  0 ; 0  0  1] ; 

YI  = [ ] ;  YY= [ ] ;  L= [ ] ;  LP= [ ] ;  LR= [ ] ;  WI= [ ] ; 
pathl=[];  pathL= [ ] ;  pathw= [ ] ;  pathW=[];  dh=[]; 

%Initial  conditions  for  the  numerical  integration  routine 
Y0=[Vwbar  gammaw  phiw  psiw  27  7  0  Twbar  Lwbar  uw  psil] ; 
for  i=l:n; 

%Build  the  solution  matrices 
YI=[YI;Y0] ; 
if  rem (nnn, 1/dt) ==0 ; 

YY= [YY;Y0] ; 
end 

%Call  the  differential  equation  solver 
[t , Y] =ode23 ( ' cf f c_sim_al ' , [0  dt] , Y0 ' ) ; 

DL= length ( t) ; 

%L  flight  disturbances  for  pitch  and  roll 
if  i<=nn; 

Vl=Vlbar-5* (1+cos (pi+ ( (i-1) *pi/ (nn-1) ))); 

gammal= (-3/ (2*57 .3) ) * (1+cos (pi+ (2* (i-1) *pi/ (nn-1) ) ) ) ;h=h- 

Vl*sin (gammal) *dt; 

phil= (-5/ (2*57.3) ) * (1+cos (pi+ (2* (i-1) *pi/ (nn-1) ) ) ) ;  Ll=Llbar/cos (phil) ; 
end 

%Create  derivatives  of  relative  position  variables 
if  i>=4; 

xv= (YI ( i , 5 ) -YI ( i-1 , 5 ) ) /dt ; 
zv= (YI (i , 7) -YI ( i-1 , 7 ) ) /dt ; 
yv= ( YI ( i , 6 ) -YI ( i-1 , 6 ) )/dt; 

ya=( ( (YI ( i , 6) -YI (i-1, 6) ) - (YI ( i-2 , 6) -YI (i-3 , 6) ) )/dtA2) ; 
end 

%Update  initial  conditions  for  next  integration  time  step 
Y0=Y(DL, : ) ; 

Y0 (8) = (Twbar/cos (Y (DL, 2) ) ) +KTP* (xc-Y (DL, 5 ) )+KTD*xv; 

YO (9) = (Lwbar/cos (Y (DL, 2) ) ) +KLP* ( zc-Y (DL, 7 ) )+KLD*zv; 

YO (10) =KpP* (yc-Y (DL, 6 ) ) +KpsiP* (Y(DL, 11) -Y(DL, 4) ) +KpD*yv+KpDD*ya ; 

%Alternate  control  law  derivatives  of  relative  position  variables 
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%if  i>=4; 

%xv=(YI(i,5)-YI(i-l,5) )/dt; 

%xa=( ( (YI (i , 5) -YI (i-1,5) ) - (YI (i-2 , 5) -YI (i-3 , 5) ) )/dtA2) ; 

%  z v= (YI(i,7)-YI(i-l,7) )/dt; 

%za= ( ( (YI ( i , 7 ) -YI (i-1,7) ) - ( YI (i-2 , 7 ) -YI ( i-3 , 7 ) ) ) /dtA2) ; 

%yv= (YI (i, 6) -YI (i-1, 6) ) /dt; 

%ya=( ( (YI ( i , 6 ) -YI (i-1, 6) ) - (YI ( i-2 , 6 ) -YI ( i-3 , 6 ) ) ) /dtA2) ; 

%end 

%Alternate  control  law  initial  conditions  for  integration 
%Y0=Y(DL, : ) ; 

%Y0 ( 8 ) = (Twbar/cos (Y (DL, 2 ) ) ) +KTP* ( zc-Y (DL, 7 ) ) +KTD*zv+KTDD*za ; 

%Y0 (9) = (Lwbar/cos (Y(DL, 2) ) ) +KLP* (xc-Y(DL, 5) ) +KLD*xv+KLDD*xa ; 

%Y0 ( 10 ) =KpP* (yc-Y (DL, 6 ) ) +KpsiP* (Y (DL, 11) -Y(DL, 4) ) +KpD*yv+KpDD*ya 
%Transform  error  signals 
%e=CWI* [xc-Y (DL, 5 ) ;yc-Y(DL,6) ;zc-Y(DL,7) ] ; 

%e=-CLI 1 *CWI* [xc-Y (DL , 5 ) ;yc-Y (DL , 6 ) ;zc-Y(DL,7) ]  ; 

%Y0=Y (DL, : ) ; 

%Y0 (8) = (Twbar/cos (Y(DL, 2) ) )+KTP*(e(l) )+KTD*xv; 

%Y0 (9) = (Lwbar/cos (Y(DL, 2) ) ) +KLP* (e (3 ) ) +KLD*zv; 

%Y0 (10) =KpP* (e (2 ) ) +KpsiP* ( Y ( DL , 11) -Y(DL, 4) ) +KpD*yv+KpDD*ya ; 

%Rotating  reference  frame  transformation 
CWIll=cos (Y0 (4) -Y0 (11) ) *cos (Y0 (2) ) ; 

CWI12=sin(Y0(4) -Y0(11) )*cos(Y0(3) )+cos(Y0(4)- 
Y0(11) ) *sin (Y0 (2 ) ) *sin (YO ( 3 ) )  ; 

CWI13=sin(Y0 (4) -Y0 (11) ) *sin(Y0 (3) ) +cos (YO (4)  - 
Y0 (11) ) *sin (YO (2) ) *cos ( YO ( 3 ) ) ; 

CWI21=-sin(Y0 (4) -YO (11) ) *cos (YO (2) ) ; 

CWl22=cos (YO (4) -YO (11) ) *cos (YO (3) ) -sin(Y0 (4)  - 
Y0 (11) ) *sin(Y0 (2) ) *sin(Y0 (3) ) ; 

CWl23=-cos(Y0(4)-Y0(ll) )*sin(Y0(3) )-sin(Y0(4)- 
Y0 (11) ) *sin (YO (2) ) *cos (YO (3) ) ; 

CWI31=sin(Y0 (2) ) ; 

CWl32  =  -cos (YO (2 ) ) *sin (YO (3 ) )  ; 

CWI33=cos (YO (2) ) *cos (YO (3) ) ; 

CWI=[CWI11  CWI12  CWI 1 3 ; CWI 2 1  CWI22  CWI23;CWI31  CWI32  CWI33]; 

%Rotating  reference  frame  transformation  (L  to  I) 

%CLIll=cos (YO (11) ) *cos (gammal) ; 

%CLI12=cos (YO (11) ) *sin (gammal) *sin (phil ) -sin ( YO ( 11 ) )*cos(phil) ; 
%CLI13=sin (YO (11) ) *sin (phil) +cos (YO (11) ) *sin (gammal) *cos (phil) ; 
%CLI21=sin(Y0 (11) ) *cos (gammal) ; 

%CLl22=cos (YO (11) ) *cos (phil) +sin(Y0 (11) ) *sin (gammal) *sin(phil) ; 
%CLI23=-cos (YO (11 ) ) * sin (phil) +sin (YO (11) ) *sin (gammal) *cos (phil) ; 
%CLI31=-sin(gammal) ; 

%CLI32=cos (gammal) *sin(phil) ; 

%CLI33=cos (gammal) *cos (phil) ; 

%CLI= [CLI11  CLI12  CLI13 ;CLI21  CLI22  CLI23;CLI31  CLI32  CLI33]; 

%Inertial  separation  and  ground  track  calculations 
if  rem(nnn, 1/dt) ==0; 

j=j+l; 

WI= [WI; [CWI* [Y (DL, 5) ; Y(DL, 6) ;Y(DL, 7) ] ]’); 

pathl= [pathl; [-Vl*sin(Y (DL, 11) )  Vl *  cos ( Y ( DL ,11))  Ll]]; 

pathL= [pathL; [sum(pathl ( : , 1) )  sum (pathl ( : , 2) )  h] ] ; 

pathw= [pathw; [-Y (DL, 1 ) *sin (Y (DL, 4) )  Y (DL, 1) *cos (Y (DL, 4) ) ] ] ; 

pathW= [pathW; [sum (pathw ( : , 1) ) +WI ( j , 1)  sum (pathw ( : , 2 ) ) -WI ( j , 2 )  h- 

WI  ( j , 3 ) ] ] ; 
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%L  flight  path  in  pitch  and  roll 
L= [ L ; Vl ] ; 

LP=  [LP;  gammal*57 . 3]  ; 

LR= [LR; phi 1  *  57 . 3 ]  ; 
end 

nnn=nnn+l; 

end 

ts= [0:1:  (n-1 ) *dt] ; 

%xyz  displacement  output 
tts=[0: .1:  (n-1) *dt] ; 

figure (1) , subplot (3, 2 , 1) ,plot (tts, YI ( : ,5) , tts,xc) ; title ( 'X 
Separation ' ) ; 

subplot (3, 2, 2) , plot ( tts , YI ( : ,8) ) ; title ( 'W  Thrust' ) ; 
subplot (3,2,3) ,plot (tts, YI ( : ,6) ,tts,yc) ;title( 'Y  Separation1 ) ; 
subplot (3, 2, 4) , plot ( tts, YI ( : ,3) *57.3) ; title ( 'W  Bank  Angle' ) ; 
subplot (3, 2,5) , plot (tts, YI ( : ,7) , tts,zc) ; title ( ' Z  Separation' ) 
subplot (3, 2, 6) , plot ( tts, YI ( : ,9) ) ; title ( 'W  Lift' ) ; 

%Velocity  change  output 

figure (2) , subplot (2 , 2 , 1 ) ,plot (ts,L( : ,1) ) ; title ( 'L  Velocity' ) ; 
subplot (2, 2, 2) , plot ( ts , YY ( : , 1) ) ; title ( 'W  Velocity' ) ; 
subplot (2, 2,3) , plot ( ts , YY ( : ,5) , ts,xc) ; title ( 'X  Separation' ) ; 
subplot (2, 2, 4) , plot ( ts , YY ( : , 8) ) ; title ( 'W  Thrust' ) ; 

%Altitude  change  output 

figure (3) , subplot (3 , 2 , 1) ,plot (ts, YY( : , 5) ,ts,xc) ; title ( 'X  Separation' ) ; 
subplot (3, 2,2) , plot ( ts , YY ( : ,8) ) ; title ( 'W  Thrust' ) ; 
subplot (3,2,3) , plot ( ts , YY ( : , 7) , ts, zc) ; title ( 'Z  Separation' ) 
subplot (3, 2, 4) , plot ( ts , YY ( : ,9) ) ; title ( 'W  Lif t ' ) ; 
subplot (3, 2,5) , plot ( ts , YY ( : ,2) *57 .3) ; title ( 'W  Pitch  Angle' ) ; 
subplot (3,2,6), plot ( ts , pathL ( : , 3 ) ) ; title ( 1  Formation  Altitude  Change ’ ) 
%Heading  change  output 

figure (4) , subplot (2 , 2 , 1) ,plot (ts, YY( : , 5) , ts,xc) ; title ( 'X  Separation' ) ; 
subplot (2, 2, 2) , plot ( ts , YY ( : ,8) ) ; title ( 'W  Thrust' ) ; 

subplot (2 , 2, 3) ,plot (ts,WI ( : , 1) ) ; title ( 'X  Separation  (Inertial) ' ) ; 
subplot (2,2,4) , plot ( ts , YY ( : , 1) ) ; title ( 'W  Velocity' ) ; 

figure (5) , subplot (2, 2 , 1) ,plot (ts, YY( : , 6) , ts,yc) ; title ( ' Y  Separation' ) ; 
subplot (2, 2,2) ,plot (ts,YY( : ,3) *57 .3) ;title( 'W  Roll  Angle'); 
subplot (2, 2, 3) ,plot(ts,WI( : ,2) ) ; title ( 'Y  Separation  (Inertial) 1 ) ; 
subplot (2,2,4) ,plot ( ts , YY ( : ,4) *57.3) ; title ( 'W  Heading  Angle'); 
figure (6) , subplot (2 , 2 , 1) , plot ( ts , YY ( : ,7) , ts, zc) ; title ( 'Z  Separation' ) ; 
subplot (2, 2, 2) ,plot (ts, YY ( : , 9) ) ; title ( 'W  Lift' ) ; 

subplot (2, 2,3) ,plot(ts, WI ( : , 3) ) ; title ( 1 Z  Separation  (Inertial) ’ ) ; 

subplot (2,2,4) , plot ( ts , YY ( : , 2 ) *57 . 3 ) ; title ( 1 W  Pitch  Angle'); 

figure (7) , subplot ( 2 , 2 , 1 ) ,plot (ts,LR( : , 1) ) ; title ( 'L  Roll  Angle' ) ; 

subplot (2,2,2) ,plot(ts,LP( : ,1) ) ; title ( 'L  Pitch  Angle'); 

subplot (2,2,3) , plot ( ts , pathl ( : , 3) ) ; title ( ' L  Lift ' ) ; 

subplot (2 , 2,4) ,plot (ts, YY( : , 11) *57 . 3) ; title ( 'L  Heading  Angle ' ) ; 

%Three  dimensional  turn  output 

figure (8) , subplot (2 , 2 , 1) ,plot3 (pathL( : ,1) ,pathL( : ,2) ,pathL( : ,3) ,pathw( 
, 1 ) , pathW ( : , 2 ) , pathW (:, 3 ),'.') ; 

title ( ' 3D  Flight  Path '); legend (' L ', 'W' ); grid  on 

subplot (2,2,2), plot (pathL ( : , 1 ) , pathL ( : , 3 ) , pathW ( : , 1 ) , pathW (:, 3 ),'.')  ; 
title (' Flight  Path  From  Behind' ); legend (' L ', 'W' ) ; 

subplot (2,2,3), plot (pathL ( : , 2 ) , pathL ( : , 3 ) , pathW ( : , 2 ) , pathW (:, 3 ),'.') ; 

title (' Flight  Path  From  The  Right  Side '); legend (' L ', 'W' ) ; 

subplot (2, 2, 4) , plot (pathL ( : , 1) ,pathL( : , 2) , pathW ( : , 1) , pathW ( : ,2) ,'.'); 
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title (' Flight  Path  From  Above '); legend (' L 'W' ) ; 

%Alternate  control  law  xyz  displacement  output 
%tts= [ 0 : .1: (n-1) *dt] ; 

% figure ( 1 ) , subplot (3,2,1) ,plot (tts, YI ( : , 5) , tts,xc) ; title ( 'X  Separation' ) 
%subplot (3,2,2) , plot ( tts , YI ( : ,9) ) ; title ( 'W  Lift' ) ; 

%subplot (3,2,3) , plot (tts, YI ( : , 6) , tts,yc) ; title ( ' Y  Separation' ) ; 

%subplot ( 3 , 2 , 4 ) , plot ( tts , YI ( : , 3 ) *57 . 3 ) ; title ( ' W  Bank  Angle'); 

%subplot (3,2,5) , plot ( tts, YI( : ,7) ,tts,zc) ; title ( 'Z  Separation' ) 

%subplot (3,2,6) ,plot(tts,YI( : ,8) ) ; title ( 'W  Thrust' ) ; 

%Alternate  control  law  Velocity  change  output 
%f igure (2 ) , subplot (3 , 2 , 1) ,plot(ts,L( : ,1) ) ;title( 'L  Velocity' ) ; 

%subplot (3,2,3) ,plot ( t s , YY (:,!)) ; title ( '  W  Velocity ' ) ; 

%subplot (3,2,4) , plot ( ts , YY ( : ,5) ,ts,xc) ; title ( 'X  Separation' ) ; 

%subplot (3,2,5) ,plot (ts, YY ( : , 9) ) ; title ( 'W  Lift ’ ) ; 

%subplot (3,2, 6) , plot ( ts , YY ( : , 8 ) ) ; title ( 'W  Thrust’ ) ; 

%  Alternate  control  law  Altitude  change  output 
%f igure (3) , subplot (3,2,1) ,plot (ts, YY( : , 5) , ts,xc) ; title ( 'X  Separation' ) ; 
%subplot (3,2,2) , plot ( ts , YY ( : , 9 ) ) ; title ( 'W  Lift' ) ; 

%subplot (3,2, 3) ,plot (ts, YY( : , 7) , ts, zc) ; title ( ' Z  Separation' ) 

%subplot (3,2,4) , plot ( ts , YY ( : , 8 ) ) ; title ( 'W  Thrust' ) ; 

%subplot (3,2,5) , plot ( ts , YY ( : ,2) *57.3) ; title ( 'W  Pitch  Angle'); 

%subplot (3,2,6) , plot ( ts , pathL ( : ,3) , ts,pathW( : , 3) , ’ . ' ) ; title ( 'Altitude 
Change ' ) ; 

%  Alternate  control  law  Heading  change  output 
%f igure (4) , subplot (2,2,1) ,plot ( ts , YY ( : , 5) , ts , xc) ; title ( 'X  Separation' ) ; 

% subplot (2,2,2) , plot ( ts , YY ( : ,9) ) ; title ( 'W  Lift' ) ; 

% subplot ( 2 , 2,3) ,plot ( ts, WI ( : , 1) ) ; title ( 'X  Separation  (Inertial) ' ) ; 
%subplot (2,2,4) , plot ( ts , YY ( : , 1 ) ) ;title( 'W  Velocity' ) ; 

% figure (5) , subplot ( 2 , 2 , 1 ) ,plot(ts,YY( :,6) ,ts,yc) ; title ( 'Y  Separation' ) ; 
%subplot ( 2 , 2 , 2 ) , plot ( ts , YY ( : , 3 ) *57 . 3 ) ; title ( ' W  Roll  Angle'); 

% subplot (2 , 2, 3) ,plot (ts, WI ( : , 2) ) ; title ( ' Y  Separation  (Inertial) ' ) ; 

% subplot (2 , 2,4) , plot ( ts, YY ( : , 4) *57 . 3) ; title ( 'W  Heading  Angle ' ) ; 

% figure (6) , subplot (2 , 2 , 1) ,plot(ts,YY( : ,7) ,ts,zc) ; title ( 'Z  Separation' ) ; 

% subplot (2,2,2) , plot ( ts , YY ( : , 8) ) ; title ( 'W  Thrust' ) ; 

%subplot (2 , 2,3) , plot (ts, WI ( : , 3 ) ) ; title ( ' Z  Separation  (Inertial) ' ) ; 
%subplot (2,2,4) , plot ( ts, YY ( ; , 2) *57 . 3 ) ; title ( ' W  Pitch  Angle'); 

%f igure (7) , subplot (2 , 2 , 1) ,plot (ts,LR( : , 1) ) ; title ( 'L  Roll  Angle' ) ; 
%subplot (2,2,2) , plot ( ts , LP ( : , 1) ) ; title ( ' L  Pitch  Angle'); 

%subplot (2 , 2 , 3 ) ,plot (ts,pathl ( : ,3) ) ; title ( 'L  Lift' ) ; 

%subplot (2,2,4) , plot ( ts , YY ( : ,11) *57.3) ; title ( 'L  Heading  Angle' ) ; 

Simulation  Equations  of  Motion 

%Dif ferenctial  Equations  of  formation  model 
function  eqn=cf f c_sim_al ( t , Y) ; 

global  mw  ml  g  rho  Cdow  Kw  Sw  Ll  phil  Yw  Vl  gammal; 
global  xc  yc  zc  KTI  KLI  Kpl; 

%  Vw=Y ( 1 ) ;  gammaw=Y ( 2 ) ;  phiw=Y(3);  psiw=Y(4);  x=Y(5);  y=Y(6); 

z=Y  ( 7  )  ; 

%  Tw=Y ( 8 ) ;  Lw=Y ( 9 ) ;  p=Y(10); 

%psil=Y (11) ; 

%Initialize  solution  vector 
eqns= [ ] ; 

%Solve  for  dynamic  pressure 
qw= . 5*rho*Y ( 1) A2 ; 
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%Rotating  reference  frame  transformation 
cl=cos(Y(2) ) *cos (gammal) *cos ( Y ( 11 ) -Y(4) )+sin(Y(2) ) *sin (gammal) ; 
c2=sin(Y (2) ) *cos (gammal) *sin(Y(3) ) *cos (Y(ll) - 
Y (4) ) +cos (gammal) *cos(Y(3) ) *sin(Y(ll)-Y(4) )  - 
cos ( Y (2 ) ) * sin (gammal) *sin(Y(3) ) ; 

c3=sin(Y(2) ) *cos (gammal) *cos (Y(3) ) *cos ( Y ( 11 ) -Y  ( 4 ) )  - 

cos (gammal) *sin(Y (3) ) *sin(Y(ll) -Y(4) ) -cos(Y(2) ) *sin (gammal) *cos (Y(3) 
%Dif ferential  equations  of  motion 
DY1= ( { Y ( 8 ) - (Cdow*qw*Sw+Kw* ( ( (Y ( 9 ) ) A2 ) / (qw*Sw) ) ) ) /mw) +g*sin(Y(2) ) ; 
DY2= ( -Y (9) *cos (Y(3) ) / (mw*Y(l) ) ) + ( (g*cos (Y (2 ) ) ) / Y ( 1 ) ) - 
(Yw*sin (Y (3 ) ) / (mw*Y (1) ) )  ; 

DY3=Y (10) ; 

DY4=( ( -Y (9 ) *sin (Y (3 ) ) +Yw*cos (Y (3 ) ) ) / (mw*Y(l) *cos (Y(2) ) ) ) ; 

DY5= ( ( -Y ( 6 ) *g*cos (Y (2 ) ) *sin(Y(3) ) - 

Y ( 7 ) *g*cos (Y (2 ) ) *cos ( Y (3 ) ) ) /Y(l) ) + ( ( Y ( 6 ) *Yw+Y ( 7 ) *Y ( 9 ) ) / (mw*Y(l) ) ) - 
Y ( 1 ) +Vl*cl ; 

DY  6  =Y ( 7 ) *Y(10)  +  ( (Y(7) *tan(Y(2) ) * (Y(9) *sin(Y(3) ) -Yw*cos (Y(3) ) ) - 
Y ( 5 ) *Yw) / (mw*Y ( 1 ) ) )+( (Y(5) *g*cos(Y(2) ) *sin (Y (3 ) ) ) /Y ( 1 ) )+Vl*c2; 

DY7=-Y ( 6 ) *Y ( 10 ) + ( ( -Y ( 6 ) *tan (Y(2) )*(Y(9)*sin(Y(3)) -Yw*cos (Y(3 ) ) ) - 
Y ( 5 ) *Y ( 9 ) ) / (mw*Y (1) ) ) +( (Y(5) *g*cos(Y(2) ) *cos(Y(3) ) ) /Y ( 1) ) +Vl*c3 ; 

%Controller  equations 
DY8=KTI* (xc-Y (5) ) ; 

DY9=KLI* ( zc-Y ( 7 ) ) ; 

DY10=KpI* (yc-Y ( 6 ) )  ; 

%Alternate  control  law  controller  equations 
%DY8=KTI* (zc-Y (7) ) ; 

%DY9=KLI* (xc-Y ( 5 ) )  ; 

%DY10=KpI* (yc-Y (6) )  ; 

%Transform  error  signals 
%e=CWI* [xc-Y (5) ;yc-Y ( 6 ) ; zc-Y (7 ) ] ; 

%e=-CLI ' *CWI* [xc-Y (5) ; yc-Y ( 6 ) ;zc-Y(7) ] ; 

%Controller  equations  for  different  turn  formations 
%DY8=KTI* (e ( 1 ) ) ; 

%DY9=KLI* (e (3 ) )  ; 

%DY1 0=KpI * (e (2 ) ) ; 

%L  heading  equation 

DY11= ( -Ll*sin (phil ) / (ml *Vl*cos (gammal ) ) ) ; 

%Solution  Vector 

eqn= [eqns  DY1  DY2  DY3  DY4  DY5  DY6  DY7  DY8  DY9  DY10  DY11] 
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